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A LEAST -SQUARE -DISTANCE CURVE-FITTING TECHNIQUE 


By John Q. Howell 
Langley Research Center 

SUMMARY 

A method is presented for fitting a function with n parameters 
y = f(Q!j,Q !25 . . •jQ'nJ to a set of N data points <(xi,yj^ in a manner that mini- 
mizes the sum of the squares of the distances from the data points to the curve. A 
differential-correction scheme is used to solve for the parameters in an iterative man- 
ner until the best fit is obtained. Two methods for finding the distances from the data 
points to the curve and a listing of the curve-fitting computer program are also given. 

INTRODUCTION 


Most of the generally used methods of fitting a curve to a set of data points mini- 
mize a function of the vertical distances from the points to the curve. For example, if 
{xijYi} is a set of N points and y = f(o!j^,Q! 2 , . . .,“m is a curve with n param- 
eters, then the method of least squares gives values of the parameters that minimize 

N 

E = ^[yi - • • -’“m Xi)]^ (1) 

i=l 

This may be done by taking partial derivatives with respect to the parameters and setting 
each of the resulting equations equal to zero; that is, 



where 

fi = • • -»“n; xi) 


0 = 1,2,. . ., n) (2) 


This set of n equations, sometimes called the normal equations, is then solved 
for the parameters. As is well known, if f is linear in the parameters, for example, a 
polynomial in x (ref. 1), a set of simultaneous linear equations merely has to be solved. 
However, in general, more complicated functions yield simultaneous equations that are 


nonlinear. In this case f may be expanded in a truncated Taylor series about a point 
in parameter space, and in this manner the nonlinear normal equations can be linearized 
and solved by iteration. The end result is a set of parameters that yields a minimum of 
equation (1). This is called either the Gauss-Newton method or the method of nonlinear 
least squares. However, when f(x) has a region where its derivative is large or when 
both Xj and y^ have similar error bounds, it may be more desirable to minimize the 
distance from each data point to its nearest point on the curve. This minimum distance 
is the same as the perpendicular distance from the data point to the curve. Scarborough 
(ref. 2) gives a method for curve fitting that minimizes the sum of the squares of these 
distances but his method is limited to first-order polynomials. Reed (ref. 3) and Kendall 
and Stuart (ref. 4) give schemes that are applicable to polynomials of higher order. 

These same methods are useful for any function in which the parameters enter in a linear 
fashion. Guest (ref, 5) describes a related technique that minimizes the perpendicular 
distance from each data point to a straight line tangent to the curve. This tangent is 
taken at the point on the curve having the same x-coordinate as the data point. 

The purpose of the present work is to derive and demonstrate the use of a curve- 
fitting technique that minimizes the least-square distances from each data point to the 
curve. The technique described herein works for a general function f and is most 
useful when the function being fitted contains regions where the slope is small as well as 
where the slope is large. It is also useful when the data points have error bounds asso- 
ciated with both the x- and y- coordinates. In the latter case, this technique implicitly 
assumes identical error for both coordinates. Generally the data points and the function 
can be scaled so that this condition is met. Other techniques that are less time con- 
suming may also be used in these situations. For example, a judicious choice of weights 
often makes possible the use of a standard least-squares procedure. However, in these 
cases a particular choice of weights seldom works for more than a few sets of data. The 
technique described herein does not have this disadvantage since it provides a fit even 
when all the weights are set equal to 1. 


SYMBOLS 


Ajij defined by equation (10) 

Bj defined by equation (9) 

Di distance from (xi,yi) to nearest point on curve, equations (4) and (16) 

Dj value of Dj using old parameters, equation (6) 
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Dj(x) distance from (xi,yi) to some point on curve, equation (17) 

E sum of squares of distances from data points to curve, equations (1) and (3) 

f(x) function to be fitted to data points; y = f(x) 

fj value of function at xj, equation (2) 

fi value of function at Xj with old parameters, equation (8) 

n number of function parameters 

N number of data points 

Wj weight associated with data point (xj,y^) 

(xi,yi) coordinates of one data point of set to which y = f(x) is being fitted 

(xi,yi) coordinates of point on curve nearest data point (xi,yi) 

q!^,q' 2 , . . .,Q!n parameters of y = f(x) 

Q!j,Q! 2 , . . .,0!n old parameters of y = f(x) during iteration to find least-square- 

distance fit 

a root-mean-square deviation of data points from curve, equation (14) 

DERIVATION AND DISCUSSION OF NEW TECHNIQUE 

There is given a set of N data points {xi,yj^ to which is to be fitted the function 
y = f^Q! 2 ,Q! 2 , . . ., 0 !^; x), where Q!i,Q! 2 j • • •jO'n are parameters. To obtain this fit the 
sum of the least-square distances (sum of the squares of the shortest distances) from 
the data points to the curve must be minimized. By using one of the techniques given in 
appendix A, the coordinates (xi,yj^^ are found for the point on the curve that is nearest 
each data point (xi,yi). Then it is desired to minimize 

N 

E = ^ wPi^ (3) 

i^l 
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where the distance from the ith data point to the curve is 

n/2 


1 

f/ - \2 ■ / 

- x2“l 

II 

Q 

(Xi - Xi) + (j 

^i-yO_ 


( 4 ) 

The weight to be associated with each point is given by wj. To minimize equation (3) it 
is necessary to solve the set of n normal equations 

N 

-^=0 = 2 y (i = l,2, . . .,n) (5) 

3 i=l ^ 


In general this is a set of nonlinear simultaneous equations. To solve the set, an 
iterative procedure sometimes called the method of differential correction can be used. 
First a Taylor series expansion of Di is made about some point (q' 1,0'2, . . in 

parameter space. Since Dj is a known function of the parameters, the expansion can 
easily be written 


Di(0!i,Q'2, 


.Q!n) = Dj 


n ^ 
V 9Di 

1 

n 

V 


4- — 

> 

^ 90!^ ^ 

2 

Lj 

k=l k 


k,Z=l 


8^Dj 

da^dai 


Ao!^. Aa^ + 


( 6 ) 


where 

Di = Di(o!j^,Q!2, . . 

and 

AQ!j^ = Q!jj - 

Now in equation (6) all higher order terms are dropped and only the terms that are linear 
in Aq! are kept. Then substitution of equation (6) into equation (5) gives 








Also, equation (4) yields 


ebj _ fj - Yj 9fj 
90-k f). 90!k 

where 

% = ■ * ’’^n’ 


0 = 1,2,. . ., n) 


( 7 ) 


( 8 ) 
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Now, by definition, 


and 


Bj = - 2 Wi(fi - 


i=l 


Ajk - 


I 

1=1 



^ A 

_fi_^ 

8Q!j 5- 8Q!j 


(9) 


( 10 ) 


and equation (7) then becomes 
n 

Bj = ^ Ajk Ao-k (i = 1, 2, . . ., n) (11) 

k=l 

The equations needed to fit the curve y = f(Q!j,Q! 2 , . . ^o the set of 

N points {xi,yi^ have now been derived. To use this procedure, a starting point in 
parameter space is first chosen and designated ^Q!j^,Q! 2 , . . .,«n)- Next the distances 
from each data point to the curve are found by using perhaps one of the techniques out- 
lined in appendix A. Then Bj and Aj^ are found from equations (9) and (10) and the 
simultaneous linear equations in equations (11) are solved for the quantities 
Aq!x, Ao! 2 , . . AQ!n. Lastly, the new parameters are obtained from 


“k= “k + Ao-k 


(k= 1, 2, . . ., n) (12) 


This set of new parameters is then used as a new starting point and the cycle repeated. 
This iteration is carried out until o!k converges (Ao!k «< Q!k) until it is obvious that 
a convergence will not be achieved. In the latter case a better starting point in param- 
eter space generally leads to convergence. It should be pointed out that the point in 
parameter space to which equations (9) to (11) converge may not be an absolute minimum, 
that is, the best fit. The end point of the process may be either a relative maximum or 
a relative minimum. The former case is rather unlikely but at any rate is easily detected 
by inspection of the value of equation (3) after each iteration. In the latter case a new 
starting point in parameter space must be chosen to see if convergence is achieved to a 
point where a smaller value of equation (3) is obtained. Unfortunately, it is in general 
difficult to tell when the best fit has been found, but once a fit sufficient for the particular 
need is located, it is not necessary to search farther. 
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EXAMPLES OF APPLICATION OF NEW TECHNIQUE 


The two examples given herein are chosen to demonstrate that for some cases 
the least-square-distance curve-fitting technique gives better results than the standard 
least-squares method. The first example arose when the author was trying to reduce 
some experimental plasma-physics data and led ultimately to the least -square-distance 
curve -fitting technique described in this paper. The second example is chosen since it is 
commonly known that standard least-squares procedures do not work well on this type 
of function. For the examples presented here the weights are set equal to one. By prop- 
erly choosing the weights it may be possible to obtain a fit with the least-squares tech- 
nique that is as good as that obtained with the least-square-distance method. However, 
for a different function and often for a different set of data points, a new set of weights 
would have to be chosen to achieve a good fit again. The least-square-distance technique 
described in this report does not have this disadvantage. 

Because of the extra computations involved in finding the closest point on the ‘curve, 
the least-square-distance method takes more computer time than the least-squares 
method. Based on the following two examples it is determined that the least-square - 
distance method is longer by a factor of approximately 2.5. 


Example I 


Example 1 is taken from the field of plasma physics where a common diagnostic 
tool is the Langmuir probe. The current versus voltage characteristic of this probe is 
given approximately by 


a 


y = Q'jx'^ + Q!2X + 0-2 + Q!4e 


5 ^ 


(13) 


where x is the voltage, y is the current, and q'j,Q' 2> • • a^re a set of adjustable 

parameters. The parameters aq and a 5 are always positive so the exponential term 
is large for x positive and small for x negative. 

When the Langmuir probe is used as a diagnostic tool, the current is typically mea- 
sured for a large set of voltage points. Then some curve-fitting technique is used to 
obtain a fit to the experimental data. The value of is of interest as the electron 

temperature can be obtained from it. This temperature is then used to calculate a par- 
ticular voltage in the region where the exponential term is small and the fitted function 
is used to obtain the corresponding current. The ion density can then be obtained from 
this current. From this description it is apparent (1) that both x and y may be in 
error and (2) the fit to the experimental data must be good both in regions where the 
exponential term is large and where it is small. 
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Figure 1 shows the fit tha 
least-squares technique. The ^ 
square deviation obtained by th« 
ure. The weighted root-mean-i 

,\'/2 


/N 









where is the vertical distance in the case of the least-squares technique and the 
shortest distance from the point to the curve for the least-square-distance method. The 
least-squares fit for x negative is not acceptable (fig. 1) but the fit is good for x posi- 
tive where the exponential term is large. 

The same set of data points is then used in a program based on the least-square- 
distance technique derived in the present paper. The result of this fit is shown in fig- 
ure 2. It is immediately apparent that the fit is much better and is in fact good enough 
to extract the desired information for the further data analysis as described earlier. 

The legends of figures 1 and 2 show that only agrees to within 20 percent. As 



Figure 2.- Least-square-distance fit used "with equation (13) and same data points 
as In figure 1. ajL = 0.05013; (X 2 ~ 0.599^; = 0.565^; = 2.398; 

05 = 2. 2611; a = 0.077- 
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would be expected from a comparison of figures 1 and 2, the coefficients of the polyno- 
mial portion of the function are in violent disagreement. The reason the least-squares 
technique does not do so well is that it finds a fit that is good in the high-slope and high- 
magnitude region at the expense of the fit in the small-slope and small-magnitude region. 

Example n 

Example II is chosen to show that the least-square -distance curve-fitting technique 
fits functions with a singularity. The function chosen is 

y = oi^x^ + Q-gX^ + Qigx + 0-4 + — ^ (15) 

In figure 3 the results of the nonlinear-least-squares curve-fitting scheme are 
shown. The values of the parameters and the weighted root-mean-square deviation of 
the points from the curve are shown in the legend of the figure. The least-square- 
distance fit of the same function to the same points is shown in figure 4. As in example I, 
the fit to the small-magnitude points is better when the least-square -distance technique 
is used while both methods give similar fits for the large -magnitude points. In 
example II the initial parameter guess for the least-square-distance method is more 
critical than usual. With a bad initial guess both distance -finding techniques described 
in appendix A sometimes achieve convergence to a point on the wrong side of the singu- 
larity. This can also happen if the fitting function has a very sharp peak, in which case 
the distance -finding scheme may achieve convergence to a point on the wrong side of the 
peak. Of course, if it is desired to fit a function of this t}rp6 to several sets of data, the 
program can be designed to alleviate this problem, but the fitting routine has to be dif- 
ferent for each particular function. 

CONCLUDING REMARKS 

In the present paper the least-square-distance curve-fitting method is derived and 
examples of its use are presented. This technique fits a function with n parameters 
y = i{oLY,oL2, - . •jQ'nj ^ ^ data points ^Xi?yi]y by minimizing the sum of 

the squares of the distances from the data points to the curve. A differential-correction 
scheme is used to solve for the parameters in an iterative manner until the best fit is 
obtained. Two examples of the use of this technique are presented, both involving func- 
tions having large slope variations. In both cases the least-squares fit is found to be 
lacking when compared to the least -square -distance fit. It is found that the least-squares 
technique fits the curve to points in the regions of large slope and large magnitude at the 
expense of the fit in regions of small slope and small magnitude. This does not happen 
for the least-square -distance method presented in this paper since the sum of the squares 
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of the distances from the data points to the curve is minimized. Hence, for functions of 
this t3^e the least-square-distance technique fits a function to a set of points more accu- 
rately than the least-squares method, unless much time is spent in customizing the 
least-squares weights to the particular function and particular set of data. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., June 10, 1971. 
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APPENDIX A 


TWO NUMERICAL METHODS FOR FINDING THE DISTANCE 
FROM A POINT TO A CURVE 


In appendix A two methods are presented for finding the distance from the data 
point to the curve y = f(x). The first method minimizes the distance from the 

curve to the data point, while the second method finds the perpendicular from the curve 
to the data point. Both these methods locate the point on the curve Jxi,f(Xi^ nearest 
the data point. The distance from the data point to the curve is then given by 


Di = 


(Yi - Yi) + (xi - xi) 




(16) 


For some very simple cases this point can be found analytically but the assumption is 
made here that f(x) is of such complexity that this is impossible. 


Method I 


The distance from some point on the curve y = f(x) 


Di(x) 




to the data point (xj^^y^) is 


(17) 


Now an x such that D^ is minimum may be found by solving 



(18) 


For the case of D^ o (for Dj^ = 0, the trivial solution is x^ = Xi and y^ = yi), it is 
seen from equation (18) that x^ must satisfy the equation 



Once Xi is found, y^^ is obtained from y^^ = f(Xi), and equation (16) is used to find Di. 
Equation (19) can be solved by any convenient method. For cases where the second 
derivative of f(x) is obtainable, the author has used the Newton-Raphson method with 
good success. It should be kept in mind that in some cases the solution of equation (19) 
may yield a Di that is a relative maximum or a relative minimum instead of the abso- 
lute minimum that is desired. Fortunately these cases are rare. 
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APPENDIX A — Concluded 


Method n 

The second method may be called the method of successive tangents and does not 
require higher derivatives of f(x). Consequently it is much more useful for complicated 
functions. To use this method a point on the curve is initially chosen near where the 
closest point is thought to be. This initial guess may be designated |xj,f(xi^ and a 
straight line fitted through this point tangent to the curve. This can be done by using 
either f(Xj) and f '(Xi) br f(xi) and f(xj + Ax^. In the latter case Ax is some 
small arbitrarily chosen increment. Once the straight line is found, a perpendicular is 
dropped to it from the data point ^Xj,yj^. A better estimate of the closest point on the 
curve is now obtained by letting the new Xj be the x- coordinate of the foot of the per- 
pendicular on the straight line. Then a second straight line tangent to the curve may be 
fitted through the new point j^^,f(Xi^. This process is repeated until two successive 
Xj's agree to within some previously chosen increment. For cases where f'(x) is 
easily obtained the author has used this scheme with good success. If the Newton- 
Raphson method is used with the first method and if f”(x) is zero these two methods 
are equivalent. 
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APPENDIX B 


COMPUTER PROGRAM FOR LEAST-SQUARE -DISTANCE TECHNIQUE 

Appendix B contains a description and listing of a least-square-distance curve- 
fitting program written in FORTRAN. The procedure for finding the distance from the 
data point to the curve is built into the curve -fitting subroutine. The method used is the 
successive-tangent method described in appendix A. The curve-fitting subroutine also 
has a d am ping procedure (ref. 6) included for increased stability. This program has 
operated satisfactorily for the author with several different functions but has not been 
tested extensively. 

Main Program 

It is felt that a description of the main program is not needed since any potential 
user has to write the main program around his own particular application. 

Least-Square-Distance Curve- Fitting Subroutine 

This subroutine assumes the existence of a linear-simultaneous-equation solver 
called SIMSOL. It is called by the statement 

Call SIMSOL(A,B,M) 

and solves the equation in M unknowns given by AX=B. The solution vector for X is 
returned in B, The curve-fitting subroutine also calls the subroutine FUNC described 
subsequently. A description of the calling procedure for the curve -fitting subroutine 
follows. 

Use: Call LSD(X,Y,W,N,AL,M,ERR,RMS) 

X'\ Vectors containing x- and y-coordinates of data points to which function 

yJ is being fitted. 

W Vector containing weight associated with each point. 

N The number of points being supplied to subroutine by main program. 

AL Vector containing values of function parameters. Initially a trial set must 

be supplied. The curve -fitting subroutine iterates and returns a better 
set. 

M The number of parameters in function being fitted. 
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APPENDIX B - Continued 


ERR An error criterion that must be supplied to subroutine by main program. 
The subroutine iterates until RMSold " RMS^ew ERR * RMS^ew 

Weighted root-mean-square deviation of data points from curve defined by 


RMS = 

where Dj^ is distance of ith point from nearest point on curve. 

(1) X,Y and W are all dimensioned 50 and hence N = 50. 

(2) AL is dimensioned 10 and hence M ^ 10. 

(3) A linear-simultaneous-equation solver must be provided as 
described previously. 

(4) A subroutine called FUNC containing information about the function 
must be supplied. An example is described next. 

Description of Subroutine FUNC 

The subroutine FUNC contains information about the function being fitted. This 
subroutine is called by the curve-fitting subroutine described previously. The subroutine 
listing included herein is used to fit equation (13) to a set of data points as shown in fig- 
ure 2 and is intended to be an example of how this subroutine may be written. 

Use: Call FUNC(X,Y,N,AL,XDER^ER) 

X Vector containing values of independent variable. 

Y Vector used to return values of dependent variable to curve -fitting sub- 

routine. For example, Y(I) must contain the value of the function 
evaluated at X(I) for I = 1 to I = N. 

N The number of X values being supplied to subroutine. If N = 1 only X(l) 

is supplied and the value of the function and its x-derivative must be 
returned in Y(l) and XDER, respectively. For other values of N, both 
Y and DER must be filled and XDER need not be calculated. 

AL Current values of function parameters being supplied to FUNC by LSD. 



RMS 


Restrictions: 
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APPENDIX B — Continued 


XDER Variable containing value of x-derivative of function evaluated at X(l). 
It need be calculated only when N = 1 . 


DER Matrix containing derivatives of function with respect to all function 
parameters, each evaluated at X(I), for I = 1 to I = N. This matrix 
must be filled by FUNC whenever N > 1. The defining equation is 


DER(K,I) = 



where is Kth function parameter. 

Restrictions: (1) X and Y are dimensioned 50 so N S 50. 

(2) AL is dimensioned 10. 

(3) DER is dimensioned (10,50). 
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APPENDIX B — Continued 


c 

c 

c 

c 

C MAIN PROGRAM 

C 
C 
C 

C THIS PROGRAM READS THE INITIAL GUESS AT THE PARAMETERS AND THE DATA 
C POINTS and then CAULS THE LEAST SQUARE DISTANCE CURVE FITTER. 

C 

C 

PROGRAM CFIT( INPUT.OUTPUT) 

D I MENS ION X { 50 ) , Y ( 50 ) , W < 50 ) ♦ AL ( I 0 ) 

ERR=1 .E-5 

1 read 1000,M, (AL( I) . I=I «M) 

1000 FORMAT! I2.X(E1 0 ) ) 

READ 1001.N.(X(I).Y<I).W<I).I*1.N) 

1001 FORMAT! I2/< 3ElO )) 

CALL LSD!X.Y«W.N, AL.M.ERR.RMS) 

GO TO 1 
END 
C 
C 
C 

c 

C least square distance curve fitting SUBROUTINE 

c 

c 

c 

c 

SUBROUTINE LSD ! X. Y . W ,N, AL , M , ERR. RMS ) 

C 

C THIS IS A LEAST SQUARE DISTANCE CURVE FITTING SUBROUTINE. 

C IT HAS BUILT IN THE SUCeSSIVE TANGENT LINE SCHEME TO FIND THE 
C DISTANCE FROM A DATA POINT TO THE CURVE. 

C IT CAULS SIMSOL(A.B.M) TO SOLVE THE LINEAR SIMULTANEOUS EQUATION 
C AX=B IN M UNKNOWNS. 

C IT CALLS FUNC TO OBTAIN INFORMATION ABOUT THE FUNCTION BEING FITTED. 
D I MENS ION X ( 50 ) . Y !50 ) , W!5o ) . AL! 1 O ) ,D ! 10.50 ) »B ( 10 ) . A ! 1 0* 10 ) 

D 1 MENS ION D I S ! 50 ) , X 1 < 50 ) « W2 ! 50 ) . I D ! 2 ) . Y 1 ( 50 > . DX ! 50 ) 

IF IPRINT=0. THERE IS NO PRINTED OUTPUT FROM THIS SUBROUTINE. 

I print=o 
IPRINT= 2 
TRMS=1 .ElO 
I D < 1 1 = 1 OHCONVERGED 
ID(2)=1 OH 
DO 2 1=1 .N 
2 DX! I )=0. 

START iteration 

DO 1 00 ITER= 1 . 1 00 

FIND CLOSEST POINTS ON CURVE 
5 DO 30 1=1 »N 
W2 ! I I=W! I ) 
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APPENDIX B - Continued 


C FIND CLOSEST POINT GIVEN AN INDIVIDUAL DATA POINT 
DO 10 K=1 .50 

IF<ABS(W( I ) l.LT.l .E-9)G0 TO 20 
XT=X{ I )+DX( I ) 

CALL FUNCIXT.Yl .1 .AL.XDER.D) 

DY=YI I )-Yl (1 ) 

DXT» f-DX( I )+DY*XDER )/< 1 .+XDER**2 ) 

C 

C THE NEXT 3 CARDS PREVENT OVERSHOOT BY DECREASING THE INCREMENT. 
C THIS GIVES INCREASED STABILITY AT THE EXPENCE OF INCREASED 
C CONVERGENCE TIME. 

TT = 2. 

IF (ABS(DX( I 1 ) .GT.l .E-4 )TT=ABS<DXTXDX( I ) ) 

DX( I )=DX( I )+DXT/< 1 .+TT*.5) 

T=ABS(OY)+ABS(DX( I ) ) 

IF (ABSIDXT ) .LT.ABS ( 1 .E-4*DX( I ) l.OR.T.LT. 1 .E-8 )GO TO 20 
10 CONTINUE 
print 1005 

1005 FORMAT!* SHORTEST DISTANCE NOT FOUND*) 

W2 ( I ) =0 . 

20 CONTINUE 

CLOSEST POINT FOUND. NOW FIND DISTANCE 
XI ( I )=X( I )+DX( I ) 

DISC I )=SQRT(DX( I )**2 + DY**2) 

IF(DIS(I).LT.l .E-8)DIS( I >=1 .E-8 
30 CONTINUE 
C 

C COMPLETE SET OF CLOSEST POINTS AND DISTANCES FOUND 
C NOW find a new SET OF PARAMETERS 
DDM=0. 

WT = 0, 

DO 40 I =1 .N 

DDM = oDM+DIS( I )*DIS( I )*W2C I ) 

40 WT=WT+W2 ( I ) 

RMS=SQRT(DDM/'WT ) 

IF ( I PR I NT. GE. 1 ) PR I NT 1 00 1 . RMS. < AL C I ) » I = 1 . M ) 

1001 FORMAT!* RMS=*G15.7/* PARAMETERS=*5G1 5. 7. ! 1 2X . 5G 1 5 .7 > ) 

44 IF ( AbS CTRMS-RMS ) .LT.eRR*RMS )GO TO 110 
TRMS=RMS 

CALL FUNCCXl .Y1 .N.AL.XDER.D) 

DO 4B K=1 .M 
DO 4=5 I =1 .N 

45 D!K. I ) = !Y! I )-Yl ! I ) )*0!K:. I )/DlS! I ) 

DO 60 K=1 .M 

DO 50 J=K.M 
A !K. J)=0. 

00 47 I =1 .N 

47 A!K.J)=A!K.J)+0!K.I)*D!J.! )*W2!I ) 

50 A!U.K)=A!K. J) 

B !K) =0. 

DO 60 1=1 .N 

60 B!K)=B!K)+D!K. I )*DIS! I )*W2! I ) 
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APPENDIX B - Continued 


C THE NEXT 6 CARDS ARE DERIVED FROM A DAMPING TECHNIQUE THAT INCREASES 
C stability as MENTIONED IN THE DESCRIPTION. IT WILL INCREASE 
C CONVERGENCE TIME TO SOME EXTENT. 

T = 0, 

DO 65 I =1 .M 
65 T=T+B< I )**3 
WW=.5*DDM/T 
DO 6T I =I fM 

67 A ( I . I )=A( I . I I+.BXWW 

THE FOLLOWING SUBROUTINE CALL SOLVES THE LINEAR SIMULTANEOUS 
EQUATION GIVEN BY AX=B WITH M UNKNOWNS. 

CALL SIMSOL(A.B*M) 

DO 70 1=1 »M 
70 AL(I )=AL(I )+B<l ) 

100 CONTINUE 

ID(1)=10HNOT CONVER 
ID(2)=1 OHGED 
110 CONTINUE 

IF{ IPRI*‘'T.GE.l )PRINT 1004, ID, RMS, ERR, ITER 
1004 FORMAT!/* SUBROUTINE LSD *2A10 /* RMS=*E15.7/* CONVERGENCE C 

*RITERION=*E15.7/* ITERATION COUNT=*I3/) 

IF( IPRINT.GE.2 )PRINT 1000, (X(n,Y{I),Xl(I),DIS(I),W<I),I=l,N) 

1000 FORMAT!* X , Y , X 1 , D I S , W=*/ < 5G1 5.7 ) ) 

RETURN 

END 


SUBROUTINE FUNC CALLED BY CURVE FITTING SUBROUTINE 


SUBROUTINE FUNC!X,Y ,N , AL , XOER , DER ) 

THIS SUBROUTINE IS CALLED BY LSD AND IS FOR THE LANGMUIR PROBE 
current VS VOLTAGE FUNCTION.' 

D I MENS ION X ! 50 ) » Y ! 50 ) , DER !10,50),AL!10) 

IF !N,NE. 1 )GO TO 50 

IF N=l, CALCULATE THE FUNCTION AND ITS X-DERIVATIVE AT X!l). 

XI =X ! 1 ) 

Y!1 )=XI*!AL!2)+AL!1 )*X1 >+AL ( 3) +AL ! 4 ) *EXP ! AL ( 5 ) *X 1 ) 

XDER =2.*AL!1 ) *X 1 +AL ! 2 ) + AL !4 ) *AL ! 5 ) *EXP ! AL ! 5 )*X 1 ) 

return 

50 CONTINUE 
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APPENDIX B - Concluded 


IF N NOT = 1* calculate THE FUNCTION AND DERIVATIVES WITH RESPECT 
TO ALL parameters AT POINTS X(I).I=1«N. 

DO loo I=1.N 
XI =X( I ) 

Tl=EXPt AL(5)*X1 ) 

T =AL(4)*T1 
X2=X1*X1 

Y( I )=AL(1 )*X2+AL<2)*X1+AL(3)+T 
DER (1«I)=X2 
DER (2» I )=X1 
DER (3»M*1« 

DER (4«I)aTl 
DER (5,I)aXl*T 
100 CONTINUE 
RETURN 
END 
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